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SUMMARY 


The  problem  of  formulating  an  orthogonal,  analytic  (i.e.,  non- 
numerical)  coordinate  system  for  use  in  a  simulation  of  the  coupled 
ionosphere/magnetosphere/solar  wind  system  which  has  been  perturbed  by  a 
high  altitude  nuclear  explosion  at  high  magnetic  latitude  has  been  inves¬ 
tigated.  Experience  with  existing  MHD  codes,  which  simulate  the  behavior 
of  mid-latitude  ionospheric  plasmas,  indicates  that  a  very  useful  coordi¬ 
nate  system  to  use  in  such  studies  is  one  which  is  aligned  with  the  geo¬ 
magnetic  field.  Unfortunately,  the  presence  of  magnetospheri c  field- 
aligned  current  systems  precludes  the  possibility  of  constructing  an  anal¬ 
ogous  orthogonal  magnetic-field-al igned  coordinate  system  suitable  for  the 
high  latitude  problem.  It  is  possible,  however,  to  characterize  an  inter¬ 
esting  class  of  possible  orthogonal  coordinate  functions  which  generalize 
the  formal  structure  of  the  dipole-field-aligned  coordinate  system  in  a 
mathematical ly  simple  way.  One  member  of  this  class  -  referred  to  above 
as  the  "zeta-coordinate  system"  -  has  been  studied  in  considerable  detail. 
This  coordinate  system  is  dipolar  in  character  at  small  distances  from  the 
origin  but  becomes  cylindrical  at  larger  distances.  It  can,  therefore,  be 
used  to  generate  computational  meshes  which  appear  to  be  better  tailored 
to  the  requirements  of  the  high  latitude  problem  than  those  which  can  be 
generated  using  standard  coordinate  systems.  The  zeta-coordinate  system 
is  double-valued.  The  nature  of  this  double-valuedness  can  be  specified 
quite  precisely,  however,  and  does  not  seem  to  offer  an  impediment  to  the 
use  of  this  coordinate  system  in  practical  applications. 
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SECTION  1 
INTRODUCTION 


For  the  past  several  years,  the  Defense  Nuclear  Agency  (DMA)  has 
hepn  pursuing  a  research  program  to  develop  a  comprehensive  understand!' ng 
of  the  phenomenology  of  high  altitude  nuclear  explosions.  The  motivation 
for  this  effort  has  been  a  recognition,  spawned  in  the  early  1960s  by  high 
altitude  nuclear  tests  (over  the  mid-Pacific  near  Johnston  Island)1,  that 
nuclear  explosions  at  ionospheric  altitudes  could  produce  widespread  and 
long  lasting  detrimental  effects  upon  radio  communication  links,  radars, 
and  optical  or  IR  sensors.  In  addition,  such  explosions  can  disrupt  the 
operation  of  electrical  equipment  through  the  phenomenon  of  electromag¬ 
netic  pulse  (EMP).  In  the  1970s,  the  primary  emphasis  of  DNA's  research 
was  directed  toward  effects  at  low  and  mid  latitudes  (low  latitudes 
because  all  of  the  high  altitudp  nuclear  tests  were  conducted  at  magnetic 
L-shells  of  about  two  or  less,  mid  latitudes  because  of  the  location  of 
CONUS).  The  scope  of  current  research  efforts  includes  examinations  of 
existing  data  from  atmospheric  nuclear  tests,  theoretical  efforts  to 
develop  new  understandings  which  go  beyond  available  data,  and  non-nuclear 
experiments . 

As  the  overall  picture  of  high  altitude  nuclear  phenomenology 
has  become  more  complete,  the  DNA  community  has  shifted  from  simply  trying 
to  sort  out  gross  effects  to  developing  a  refined  picture  of  high  altitude 
nuclear  explosions  and  the  manner  in  which  they  interact  with  their  sur¬ 
rounding  environment.  With  this  evolution  in  thinking  has  come  a  recogni¬ 
tion  that  the  high  altitude  environment  at  high  magnetic  latitudes  (the 
auroral  oval  and  polar  cap  regions)  is,  in  a  number  of  respects,  quite 
different  from  that  at  mid  and  low  latitudes. 

1  enevious  rAoe 
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From  the  nuclear  weapons  effects  point  of  view,  the  uniqueness 
of  the  high  latitude  ambient  environment  stems  from  the  departure,  at 
polar  latitudes,  of  the  earth's  magnetic  field  from  a  dipolar  geometry, 
(Refer  to  Figure  1.)  The  key  issues  are  the  physical  processes  which  are 
linked  to  the  highly  distorted  geomagnetic  field.  At  low  and  mid  magnetic 
latitudes  (L-shells  of  -12  or  less),  the  field  lines  are  closed  and 
fairly  closely  approximate  dipole  field  lines.  These  field  lines  lie 
within  the  plasmapause  and  are  shielded  by  the  magnetosphere  from  direct 
exposure  to  the  solar  wind.  In  contrast,  magnetic  field  lines  originating 
in  the  polar  regions  (high  L-shells)  extend  out  into  (and  perhaps  through) 
the  magnetosphere  where  they  are  exposed  to  the  influences  of  the  solar 
wind.  At  high  latitudes,  energy  delivered  to  the  magnetosphere  by  the 
solar  wind  can  be  transferred  to  the  polar  ionosphere  by  electrical  cur¬ 
rents  which  descend  along  polar  magnetic  field  lines  from  the  magneto¬ 
sheath.  These  are  the  Birkeland  currents.2  Studies3  by  specialists  in  the 
theory  of  current  driven  plasma  instabilities  suggest  that  these  currents, 
which  are  unique  to  high  latitude  field  lines,  may  drive  plasma  instabil¬ 
ities  which  lead  to  ionospheric  structure  at  scale  sizes  which  can  affect 
radio  and  radar  transmissions  over  a  broad  frequency  spectrum.  This  iono¬ 
spheric  plasma  structure  may  bear  similarity  to  field-aligned  structure 
which  develops  in  plasmas  produced  by  high  altitude  nuclear  explosions. 

In  recognition  of  the  significance  of  such  naturally  occurring 
effects,  and  in  view  of  the  direct  connection  of  such  effects  to  the  high 
altitude  nuclear  phenomenology  problem,  DNA  has  undertaken  a  research 
program  to  acquire  in-situ  data  at  polar  latitudes.  The  program  has 
included  the  WIDEBAND  Satellite  Program1*  and,  more  recently,  the  HILAT 
Satellite  Program.5  Each  of  these  programs  involved/involves  active 
ionospheric  probes  from  polar  orbiting  spacecraft.  DNA  has  also  pursued 
theoretical  efforts  to  understand  the  basic  phenomenology  leading  to 
disturbances  in  the  ambient  ionospheric  environment  and  to  understand  the 
connections  c*  the  physics  of  that  phenomenology  to  the  nuclear  effects 
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problem.  Data  from  the  WIDEBAND  satellite  and  the  associated  theoretical 
interpretations  of  that  data  have  had  a  direct  influence  on  the  current 
state  of  understanding  of  propagation  effects  resulting  from  structured 
plasma  environments.6  A  similar  wealth  of  information  from  the  HILAT 
satellite  is  anticipated. 

There  are  other  reasons  to  he  interested  in  high  latitude 
phenomenology.  Note  that  the  solar  wind  energy  flux  incident  on  the 
magnetosphere  is  equivalent  to  ~10  KT/hour  to  '10  MT/hour,  depending 
on  the  state  of  magnetic  activity.  Under  naturally  occurring  conditions, 
only  a  small  fraction  of  this  energy  flux  is  directly  coupled  into  the 
magnetosphere  system.  However,  it  has  been  suggested  that  a  high  altitude 
nuclear  explosion  at  polar  latitudes  could  possibly  alter  the  ambient 
magnetospheri c  current  patterns  in  such  a  way  as  to  deliver  over  a  hroad 
area  of  the  polar  ionosphere  an  energy  input  equivalent  to  perhaps  mega¬ 
tons  per  hour  of  additional  nuclear  explosions.7  The  sources  of  this 
energy  would  be  twofold:  i)  energy  stored  in  the  currpnt  systems  within 
the  magnetosphere,  am  ii)  energy  delivered  to  the  magnetosphere  by  the 
solar  wind  and  directly  transmitted  to  the  polar  ionosphere.  The  proposed 
initiator  of  this  process  and  conduit  for  the  energy  flux  is  the  nuclear 
p.ume,  a  magneti c-f ield-al i gned  column  of  high  density  plasma  created  by  a 
high  altitude  nuclear  explosion,  which  may  reach  several  tens  of  thousands 
of  kilometers  from  the  polar  ionosphere  up  into  the  magnetospheric  cur¬ 
rents  . 

At  present,  this  concept  has  not  been  explored  by  theoretical  or 
experimental  means.  The  equations  which  are  thought  to  describe  the 
coupled  ionospbere/i  iagnDtosphere/solar  wind  system  under  ambient  condi¬ 
tions  are  complex,  and  the  magnetic  field  geometry  in  which  these  equa¬ 
tions  need  to  be  rolved  is  formidable.  In  addition,  the  extra  complica¬ 
tions  introduced  by  a  high  altitude  nuclear  explosion  render  the  nuclear 
effects  problem  euite  difficult. 


A  first  step  in  attacking  this  problem  is  to  develop  a  suitable 
framework  in  which  to  perform  calculations.  Due  to  the  complexity  of  the 
problem,  it  is  assumed  that  the  theoretical  approach  will  involve  numeri¬ 
cal  calculations.  This  report  describes  a  specialized,  three  dimensional, 
orthogonal  coordinate  system  which  has  been  developed  for  the  purpose  of 
numerically  modeling  the  earth's  magnetosphere  and  the  interaction  of  high 
altitude  nuclear  explosions  with  magnetospheric  currents.  Because  the 
magnetosphere  spans  an  incredible  volume  by  terrestrial  standards,  and 
because  numerical  simulations  of  it  must  treat,  with  some  sensible  degree 
of  spatial  resolution,  features  with  a  wide  range  of  characteristic  scale 
sizes,  the  well  known  coordinate  systems  (Cartesian,  cylindrical,  and 
spherical)  are  easily  shown  to  present  formidable  problems.  If  one  wants 
to  resolve  simultaneously  "fine  scale"  features  near  the  earth  (e.g., 
Rirkeland  currents  and  associated  ionospheric  processes),  some  measure  of 
detail  in  the  bow  shock  and  magnetopause  regions  (e.g.,  coupling  mechan¬ 
isms  by  which  solar  wind  energy  is  transferred  to  the  magnetospheric  cur¬ 
rent  systems),  and  very  large  gradient  length  features  associated  with  the 
tail  region  (which  stretches  many  tens  of  earth  radii  in  the  antisunward 
direction),  then  a  coordinate  system/computational  grid  combination  with 
enough  flexibility  to  distribute  grid  points  or  cells  in  an  appropriate 
fashion  with  economy  of  computer  resources  is  required. 

When  one  sets  out  to  design  a  computational  grid  within  which  to 
model  the  magnetosphere,  one  rapidly  finds  that  direct  application  of  the 
well  known  coordinate  systems  leads  to  an  enormous  number  of  grid  cells 
which  far  exceeds  the  high  speed  memory  capacity  of  any  present-day  com¬ 
puter.  Furthermore,  none  of  these  coordinate  systems  naturally  matches 
all  of  the  boundary  surfaces  or  important  features  of  the  magnetosphere. 
For  example,  the  spherical  coordinate  system  is  ideal  for  modeling  the 
ionosphere  and  features  near  the  earth  (provided  one  makes  it  a  geocentric 
coordinate  system),  but  far  from  the  earth,  spherical  coordinate  surfaces 
don't  match  the  shape  of  the  outer  portions  of  the  magnetosphere  or  con- 
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without  loss  of  generality.  This  is  the  required  conclusion.  To  estab¬ 
lish  the  second  part  of  Theorem  2,  suppose  R ''  =0  so  that  Rj  =  r+a  j,  modulo 
an  irrelevant  multiplicative  constant.  The  analysis  above  regains  valid 
up  to  and  including  Equation  (3-16),  hut  now  it  is  possible  that  k3  does 
not  vanish.  Assume  this  happens.  Then  none  of  the  k.'s  can  vanish. 

Using  Equations  (3-9)  and  (3-10)  one  obtains 


and 


R  '  r  +  , 

r2-  -  h  (-7^) 


R  1 

_1,  *  k. 
P, 


C-V-1) 


( 3 -27a ) 


(3-?7h) 


Inserting  these  expressions  into  Equation  (3-12)  results  in 


k  i  k  2 


f  r  + 


a-l)2  =  >- 


(3 -2b) 


Since  t.Ms  true  for  all  r,  aj  must  vanish.  To  determine  <t>2  and  $3,  use 
Equation  (3-13): 


or 
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...1 


<t> 

_2 

4>‘ 


-k4  /  d$ 

<t> 

e  2 


(3 


(3-30) 


Mote  that  Equation  (3-13)  is  the  only  constraint  upon  the  functions  0.  so 
that  is  arbitrary  here.  Without  loss  of  generality,  then,  one  may 
redefine  t>2  ♦  '.  This  leads  to 
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kl  R1 
r2  R  | 


(3-21) 


or 


r  Ri dr 
r2  R' 


up  to  an  irrelevant,  nul  tipi  i  cat  i  ve  constant.  Also 


O' 


(3-2?) 


or 


J-  =  k  l  cot  8 


q2  =  sin  le 


(3-23) 


(3-2A) 


By  Equation  (3-11),  R3  and  o3  are  both  constant.  To  determine  $2  and  $3 
consider  Equation  ( 3-1 3h ) .  One  of  $2  and  <d  must  vanish.  If  vanishes 
then  the  c's  no  longer  define  a  three  dimensional  coordinate  system,  so 
that,  in  fact  4/  must,  vanish.  Thus,  $2  must  he  constant,  leaving  4>3  arbi¬ 
trary.  To  summarize,  it  has  been  shown  that 


Rj_  dr 

J  r-2 


R' 


A2 


=  e 


l  sin  e] 


lki 


As  =  M<t>) 


(3-25a) 


( 3-25b ) 


up  to  irrelevant  constants.  Again  arguments  given  in  Section  2  imply  that 
A2  and  A 3  may  he  defined  by 


j  h  dr 


2  R' 


A2 


sin  e 


( 3-26a ) 


and 


A3 


(3-26h) 
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Thp  requirement  R"  *  0  may  now  be  used  t.o  show  that  k  3  =  0.  From  (3-9h) 
and  (3-10h) 


R2R3  =k,  k,  r*i  ]2 

kl  k2  W;]  * 


(3-17) 


Substituting  into  (3-12b),  there  results 


kxA 2  (?I)2  =  k3 
r2  *3 


(3-18) 


/k,  k,  R ,=  /k,  r  R‘ 
12  1  3  1 


(3-19) 


Differentiating  both  sides  and  using  (3-16),  this  becomes 


/k R  ‘  =  A  ,  R  *  +  A  ,  rR " 
3  1  3  1  3  1 


(3-20) 


which  implies  k,  =  0  since  R"  *  0.  Of  course,  (3-16)  new  also  implies 
that  k3,  k4  and  k^2  all  vanish.  From  (3-9),  (3-10),  and  (3-11)  kj  and  k2 
cannot  both  vanish,  however,  for  then  {R ,  O.1 ,  }  would  vanish  for  one  of 

i  =  2  or  3.  This  would  mean  that  q.  =  constant  for  i  =  2  or  3,  and  {q. } 
would  no  longer  define  a  three  dimensional  coordinate  system.  Thus,  there 
are  two  alternati ves . 


Either  k2  =  k3  =  k4  =  0  with  kt  *  0  or 


k3  =  k4  =  0withk2*0. 


The  symmetry  of  form  present  in  Equations  (3-9)  and  (3-10)  implies  that 
the  coordinate  systems  determined  by  each  of  these  alternatives  are  the 
same  -  one  may  consider  only  the  first  alternative  without  loss  of 
generality.  From  Equation  (3-9)  then 
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and 
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Equation  (3-11)  in  turn  implies 
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The  right  hand  side  of  Equation  (3-12)  may  also  he  separated 


sin2e 


2 _ 3  +  k. 

Q2  °3 


_2 _ a 
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(3-1 3a ) 


=  k, 


( 3 - 1 3h ) 


All  the  k.  above  are  at  this  stage  arbitrary  constants.  The  proof  is 
completed  by  showing  that  some  of  these  constants  must  vanish  and  by  using 
this  information  to  delimit  t'  e  form  of  the  various  derivatives  occurring 
above.  Multiplying  (3-9)  by  (3-10)  one  obtains 


^a 

°2  °3 


=  k  j  k  2  cot  2e  . 


(3-14) 


Substituting  this  expression  into  (3-13a)  results  in 


sin2e  f k  j  k  2  cot  2o  +  k  3]  =  k4  . 


(3-15) 


The  validity  of  Equation  (3-15)  for  arbitrary  9  implies  the  following 
important  relationship  among  the  separation  constants 


^3=^4=^1^2 


(3-16) 


2U 


I  I  ■' 


Vr  •  7 r,,  =  K  1  P  '  CCS  9  0,  <t, 
s  »  *  3  ;  -j  id 


-7  R  1»  3  *  0;  *3 

r  7 


VC2  *  ^3  =  R'  R'Cc,  °J 


+  °’,  *2  *> 


7  7  R;  R  3  n,  *'  *’ 

r7  sin70  2  '■  < 


Dividinc  both  sides  of  (3-6)  hy  {  -1-^2  cos  9  02  $2I  flnd  noting  that  if  a 
function  of  r  alone  is  equal  to  a  function  of  fi  alone  for  all  r  and  e, 
then  these  two  functions  must  both  be  equal  to  a  constant,  one  obtains 


r*  R 1  R1 
_ 1 _ 2 

*1  R2 


-  tan  9  -i 

0, 


( 3-9a ) 


Similarly,  (3-7)  and  (3-8)  imply  respectively 


(3-9b) 


r2  r;  r;  0l 

- 1 — 3  =  tan  e  — 3 

R 1  R  3  Q3 


(3-10a) 
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(3-10b) 


I 


if  R"  *  0.  (Here  the  primes  imply  differentiation.)  If  R'j  =  0,  then  the 
following  set  of  coordinate  functions  is  the  only  additional  possibility: 


=  r  cos  0  , 


(3-4a) 


C2  =  r  sin  9  *2U)  . 


<f2d<)) 

-J 

53  =  r  sin  9  e  2 


Here  <t2($)  is  an  arbitrary  function  such  that  V  *  0. 


Proof:  The  proof  is  accomplished  by  using  the  orthogonality  requirement 
to  develop  a  system  of  simultaneous  partial  differential  equations  for  the 
R . ,  0^  and  $ . .  These  equations  may  then  be  solved  in  terms  of  a  set  of 
separation  constants  k^  to  obtain  the  required  result.  Proceeding  then, 
the  orthogonality  requirement 


Vci  •  vCj  =  0  , 


1  *  J  » 


implies 


•  Vc2  =  R[  cos  9  02  $2 


R l R 2  sin  9  O'  *2 
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SECTION  3 

A  CLASS  OF  POSSIBLE  COORDINATE  FUNCTIONS 


In  thi?  section,  t wo  results  will  he  presented  regarding  the 
possibility  of  formally  generalizing  the  mathematical  structure  of  the 
dipole-f ield -al i gned  coordinate  system  described  in  Section  2.  Many 
possible  generalizations  might  he  explored.  However,  for  the  sake  of 
analytic  t ractahi 1 ity,  attention  here  will  be  restricted  to  coordinate 
functions  of  the  form 

c  =  P(r)  0(e)  «U)  (3-1) 

where  r,  9,  <t>  are  the  usual  spherical  coordinates.  Coordinate  functions 
of  this  form  will  he  referred  to  as  "separable  coordinate  functions."  We 
present  two  important  theorems. 

Theorem  ? 


The  most  general  set  of  "separable  coordinate  functions," 
including  the  coordinate  function 

?!  =  R  ylr)  cos  e  (3-2) 

and  satisfying  the  requirements  of  orthogonality,  is  given  by 


?!  =  R 1 (r)  cos  0  , 

( 3-3a ) 

r  2P  ' 

?2  =  e  1  sin  6  , 

( 3  -  3  h ) 

<3  "  4 
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(  3  -  3  c ) 


Thus,  a  magnetic  field  R,  obeying  the  usual  Maxwell  equations,  cannot 
possess  a  pseudo-potential  unless  B  •  Jtotai  vanishes,  where  J  .j  is  the 
total  current  associated  with  B.  Therefore,  the  field-aligned  currents  in 
the  magnetosphere  prohibit  the  existence  of  a  field-aligned  coordinate 
system. 


Even  though  the  presence  of  field  aligned  currents  forbids  the 

existence  of  an  orthogonal  coordinate  system  aligned  with  the  magneto- 
♦ 

spheric  R  field,  it  may  be  still  possible  to  design  a  coordinate  system 
which  is  more  faithful  to  the  particulars  of  the  earth-solar  wind  geometry 
than  one  of  the  standard  coordinate  systems.  Such  a  coordinate  system 
might,  for  example  possess  a  dipolar  or  spherical  character  near  the  sur¬ 
face  of  the  earth  but  transition  to  a  more  cylindrical  behavior  at  large 
distances  from  the  earth  in  order  to  conform  to  the  elongated  configura¬ 
tion  of  the  magnetosphere  (c.f.  Figure  1).  The  remainder  of  this  report 
is  dedicated  to  an  exploration  of  this  idea. 
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coordinate  systpn  defined  by  Equations  (2-4)  permits  a  ready  decomposition 
of  the  plasma  motion  into  components  which  arp  either  parallel  or  trans¬ 
verse  to  this  magnetic  field.  Fxperipnce  has  shown  that  such  a  decompo¬ 
sition  is  crucial  in  simulating  the  detailed  aspects  of  plasma  behavior 
relevant  to  high  altitude  nuclear  phenomenology  while  meeting  the  con¬ 
straints  of  available  computer  resources. 

Given  this  experience,  it  seems  appropriate  to  begin  attacking 
the  problem  with  which  this  report  is  concerned  by  asking  whether  or  not 
it  is  possible  to  construct  an  orthogonal  coordinate  system  aligned  with 
the  non-dipolar,  solar-wind-distorted  magnetic  field  of  thp  earth's  mag¬ 
netosphere  (c.f.  Figure  1).  Interestingly  enough,  it  can  be  demonstrated 
that  the  presence  of  field-aligned  currents  within  the  magnetosphere 
rigorously  precludes  the  possibility  of  such  a  construction.  The  validity 
of  this  statement  is  directly  implied  by  the  following  necessary  condition 
for  the  vector  field  F(r)  to  possess  a  pseudo-potpntial .  9 

Theorem  1 


If  the  vector  field  F ( r )  has  a  pseudo-pot enti al  Y(r),  then 

-►  -► 

F  •  VxF  must  vanish. 

Proof:  F  •  VxF 

* 

=  u(r)  VY  .  Vxg(r)VY  (2-6) 

4-  -► 

=  u(r)  VY  •  [Vg(r)  x  VY  +  y(r)  VxVY]  (2-7) 

=  0 
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for  some  i.  If  Fquation  (2-3)  holds,  u(0  I?  said  to  ho  an  integrating 
factor  and  c..  a  pseudo-potent i a  1  for  the  field  F. 

An  interesting  examplp  of  these  ideas,  important.  for  MHD  Simula 
tions  of  mid-latitudp  ionosphpric  plasmas,  is  thp  geomagnet ic-dipole- 
field-al ignpd  roordinatp  system.  This  coordinatp  system  is  defined  hy  tre 
coordinate  functions*  a,  B,  Y: 

cos  a  =  cos  9  ,  (2-4a) 

R 

sin  B  =  (— )1/2  sin  e  ,  (?-4b) 

R 

and 

Y  -  *  •  (2-4c) 

Here,  Rp  is  the  radius  of  rhp  earth  and  (R,  9,  $)  are  geocentric  spherica’ 
coordinates  with  the  colatitude  9  being  measured  rrom  the  geomagnetic 
pole,  a,  B,  and  7  satisfy  the  orthogonality  conditions  (2-2).  In  addi¬ 
tion,  it  may  he  easily  verified  that  the  gradient  of  a  is  proportional  to 
the  geomagnetic  dipole  field,  ^jp0]p*  given  hy 

P  ..  ,  =  -.31  (— l3  [2  cos  9  r  +  sin  e  el  (gauss)  (2-b) 

dipole  p  J  1  1 

where  r  and  8  are  radial  and  colatitudinal  unit  vectors.  For  an  iono¬ 
spheric  plasma  moving  in  the  dipole  magnetic  field  of  the  earth,  the 


Notice  that  thp  sets  of  coordinate  functions  (£.  }  and  (f.  (l)|  (no  sum 
on  i)  describe  equivalent  coordinate  systems  since  the  gradient  of  r,. 
and  that  of  f.(£.)  are  parallel.  This  freedom  of  definition  has  been 
used  in  Equation^  (2.4)  to  define  coordinate  functions  a,  B,  Y  which 
are  measured  In  radians. 
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SECTION  2 

SOME  PRELIMINARY  CONCEPTS 


As  a  prelude  to  a  detailed  discussion  of  the  problem  with  which 
this  report  is  concerned,  it  may  he  useful  to  review  a  few  basic  defini¬ 
tions  regarding  the  concept  of  a  coordinate  system.  Following  Reference 
9,  a  coordinate  system  is  a  threefold  family  of  surfaces  with  defining 
equations 


(r )  =  (constant )i 


1  <  i  <  3  , 


which  may  be  inverted  to  yield  r  as  a  function  of  the  q  's.  The  's  are 
referred  to  as  coordinate  functions  and  the  intersection  of  two  of  the 
surfaces  defined  in  (2-1)  is  called  a  coordinate  line.  The  coordinate 
system  defined  by  {^.  }  (=  {ci.C2»C3})  is  said  to  be  orthogonal  if 


Vci  •  =  0  ,  i  *  j  .  (2-2) 

In  practice,  the  property  of  orthogonality  is  an  important  simplifying 
feature  of  a  coordinate  system.  For  this  reason,  as  stated  in  the  intro¬ 
duction,  only  orthogonal  coordinate  systems  will  be  considered  in  this 
report.  One  last  concept  which  should  he  mentioned  at  this  point  is  that 
of  a  field-aligned  coordinate  system.  A  coordinate  system  {5.  ]  is  said 

to  be  aligned  with  a  vector  field  F(r)  if  Vc.  is  parallel  to  F  for  some 

*  ♦ 

i.  This  means  that  there  must  exist  a  scalar  function  p(r)  such  that 


u(r)  Vci  =  F (r ) 
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Given  that  the  unusual  geometry  of  the  magnetosphere  is  a  cause 
for  difficulty,  one  might  he  tempted  by  the  idea  of  a  magnetic-field- 
aligned  coordinate  system.  This  approach  will  he  explained  later  in  this 
report.  The  idea  is  quite  attractive,  but  as  will  be  shown,  the  differen¬ 
tial  geometry  of  such  an  idea  is  incompatible  with  the  physics  of  the 
magnetosphere. 

The  coordinate  system  developed  in  this  report,  the  "zeta  coor¬ 
dinate  system",  has  been  designed  to  have  desirable  geometric  properties 
in  the  near-earth  region  (i.e.,  it  can  match  selected  features  in  and 
above  the  ionosphere)  and  to  transition  gracefully  to  a  coordinate  system 
which  allows  for  simple  exterior  boundary  surfaces.  The  coordinate  sur¬ 
faces  of  this  system  can  be  made  to  conform  closely  to  dipolar  surfaces 
near  the  earth  (so  they  look  like  a  field-aligned  system  there).  These 
coordinates  retain  enough  flexibility  to  permit  the  user  to  orient  the 
dipole  axis  arbitrarily  relative  to  the  earth-sun  axis.  Far  from  the 
earth,  the  coordinate  system  approaches  a  cylindrical  system  with  the  axis 
along  the  earth-sun  direction.  The  zeta  coordinate  system  is  orthogonal. 

The  following  sections  of  this  report  explain  the  theoretical 
basis  for  this  coordinate  system.  In  Section  2,  a  brief  review  of  some 
basic  mathematical  concepts  relevant  to  the  study  of  coordinate  systems  is 
presented.  In  Section  3,  an  attempt  is  made  to  generalize  the  mathemati¬ 
cal  form  of  the  dipole-field-aligned  coordinate  system8  which  has  been 
used  previously  for  modeling  low  and  mid  latitude  nuclear  effects.  This 
results  in  a  simple  characterization  of  an  interesting  class  of  possible 
coordinate  functions.  Next,  in  Section  4,  one  member  of  this  class,  the 
zeta  coordinate  system,  is  examined.  This  system  is  well  suited  to  per¬ 
forming  simulations  of  magnetospheric  physicr.  Finally,  the  results  are 
summarized  in  Section  5. 
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farther  in  the  antisunward  direction  than  in  the  sunward  direction  without 
any  problem.  However,  near  the  earth  the  coordinate  surfaces  do  not  match 
the  geometry  of  the  ionosphere.  In  order  to  resolve  details  in  the  iono¬ 
sphere,  the  cylindrical  cells  must  be  relatively  small  (-50  km  on  a 
side,  for  example).  Direct  extension  of  the  resulting  cylindrical  coordi¬ 
nate  surfaces  to  magnetospheric  distances  leads  to  a  large  number  of 
computational  cells  which  are  much  smaller  than  is  appropriate.  This 
means  that  one  must  resort  to  numerical  "fixes"  to  try  to  make  the  simula¬ 
tion  physically  sensible  in  an  important  portion  of  the  problem  (the  iono¬ 
sphere),  subject  to  the  constraint  of  a  limited  computer  memory  (i.e.,  a 
limited  number  of  grid  cells). 

A  third  candidate  is  Cartesian  coordinates.  The  situation  is 
similar  to  that  of  cylindrical  coordinates.  Far  from  the  earth,  the 
boundary  conditions  are  easily  implemented,  but  near  the  earth,  coordinate 
surfaces  do  not  match  the  natural  geometry  of  the  problem.  Closely  spaced 
coordinate  surfaces  near  the  earth  translate  into  excessive  numbers  of 
unnecessarily  small  grid  cells  well  out  into  the  magnetosphere. 

It  is  worth  noting  that  so  far  only  orthogonal  coordinate  sys¬ 
tems  have  been  considered.  This  has  been  deliberate.  In  order  to  insure 
that  the  differential  or  integral  equations  that  are  to  be  solved  to  simu¬ 
late  the  magnetosphere  remain  tractable,  we  have  chosen  to  require  ortho¬ 
gonality  of  any  candidate  coordinate  system.  Experience  has  shown  that 
the  effort  required  to  implement  on  a  computer  complex  equations  in  a  non- 
orthogonal  coordinate  system  can  become  unreasonably  large.  In  addition, 
we  impose  the  requirement  that  the  candidate  coordinate  system  be  analyti¬ 
cally  generated.  This  requirement  allows  one  to  perform  a  fair  amount  of 
exploration  of  the  physics  equations  before  going  to  the  computer. 
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form  to  the  streamlines  of  the  solar  wind.  If  one  wants  the  hounding  sur¬ 
faces  of  the  computational  space  to  be  simple  coordinate  surfaces  (spheri¬ 
cal  in  this  case),  then  a  grid  which  extends  well  out  into  the  nagnetotail 
(say  40  earth  radii)  will  also  extend  well  out  (40  earth  radii  to  he 
exact)  in  the  sunward  direction.  Unf ortunately ,  that  means  slightly  less 
than  half  of  the  entire  grid  volume  will  lie  sunward  of  the  bow  shock  in 
the  zone  of  unperturbed  solar  wind.  For  many  problems,  this  situation 
represents  a  tremendous  waste  of  computer  resources  (storage  and  central 
processor  time).  This  problem  can  be  overcome  by  programming  the  computer 
to  chop  out  or  ignore  grid  cells  in  uninteresting  regions,  but  only  at  the 
expense  of  computer  code  simplicity. 

Also  note  that  if  one  wants  spatial  resolution  of  2  earth  radii, 
for  example,  at  a  distance  of  20  earth  radii,  the  angular  separation  of 
radial  grid  lines  needs  to  he  0.1  radians.  Therefore,  in  a  simple-minded 
application  of  spherical  geometry,  these  radial  lines  define  cells  in  the 
ionosphere  which  have  dimensions  on  the  order  of  bhO  km  —  far  too  large 
to  be  useful.  Conversely,  choosing  the  cell  dimensions  and  radial  coordi¬ 
nate  surface  spacing  according  to  ionospheric  criteria  leads  to  numerous 
cells  at  large  distances  from  the  earth  which  are  inappropri ately  small. 

These  considerations  lead  to  the  conclusion  that  spherical 
coordinates  are  not  particularly  well  suited  to  the  magnetospheric 
problem.  Numerical  calculations  in  a  spherical  coordinate  system  would 
require  substantial  effort  just  to  define  an  acceptable  computational 
grid. 


It  is  appropriate  to  next  investigate  cylindrical 
Assume  the  cylindircal  axis  lies  along  the  earth-sun  line, 
may  wish  to  convince  himself  that  other  orientations  are  of 
ity,  at  best.)  Then  boundary  conditions  far  from  the  earth 
magnetosphere)  are  simple  to  implement,  and  the  grid  can  be 


coordinates . 
(The  reader 
1 imited  uti  1  - 
(outside  the 
extended  much 
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*3 


/  d<|, 


•<1  *2 


( 3-3 1  a ) 


-k2 


=  e 


( 3- 3 1  b ) 


Summarizing,  the  following  expression  for  the  5.  are  implied  by  the  equa¬ 
tions  above: 


Cj  =  r  cos  6  , 


C2  =  (r  sin  9  $2) 


and 


?3  =  (r  sin  0 


which  is  equivalent  to  the  desired  result. 


Theorem  3 


( 3- 32a ) 

( 3-32b ) 


(3-32c) 


The  only  set  of  coordinate  functions  {c.  }  satisfying  orthogonal¬ 
ity  and  separability  and  such  that 

Cj  =  Rj(r)  sin  0  cos  $  (3-33) 

is  the  set  of  cartesian  coordinate  functions. 

Proof:  The  orthogonality  and  separability  requirements  may  be  used  as  in 
Theorem  2  to  develop  a  system  of  simultaneous  partial  differential  equa¬ 
tions.  These  equations  are  then  solved  in  terms  of  a  set  of  separation 
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From  (3-34)  and  (3-35) 


r2  R'  R!  o! 


l  i  “i  ,  $! 

+  o'  cot  9  -  ~  tan 

1  i  j  sin  9 

i 


(3-37) 


for  i  =  2,  3  .  This 


implies 


r2  R  •  R  ' 

_ L__i  =  k 

R,  R  i 
1  i 


(3-38) 


r  cot  0  -  _1  =  .  k 

ui  sm2  e  <t.  Ki  • 


(3-39) 


E<w'>»  (3-39)  l"ay  he  separated  to  prodpce 


Lr  cos  0  +  ki  si"  0]  sin  e  =  k 

1  C+I 


(3-40) 


iT  tan  *  =  +  k2+T  * 


turning  to  Equation  (3-36)  one  finds 


(3-41) 


R2  r3  02  03  sTn^e^tT" 


(3-42) 
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This  implies 


r2  R 1  R' 

- 2— a  =  k6 

r2  r3 


(3-43) 


and 


$  <t 

2  3 


°2  S3 


©2  03  sinz9  «2  <&3 


(3-44) 


The  last  equation  may  he  further  separated  to  produce 


,0‘  O' 

-2—1  +  k 

02  03 


6J  Sin  0 


2«  = 


(3-45) 


and 


S' 

-2—1  =  +  k7  . 
*2  °3 


(3-46) 


The  above  equations  can  now  he  solved  for  the  various  k.  and  the  required 
result  obtained.  First,  note  that  Equation  (3-41)  implies 


-2—1  tan 2<t>  =  k5  k„ 
*2  *3 


(3-47) 


With  the  help  of  Equation  (3-46)  this  becomes 


k?  tan24,  =  k5  k4 


(3-48) 


which  implies 


$'  $'  4>‘ 

Thus,  hy  Equation  3-46,  one  of  — 2  ,  —• 1  must  vanish.  Suppose  —2  =  0. 

*2  *7,  *2 

Then,  hy  Equations  (3-39)  and  (3-38) 


—  =  -  k2  tan  9 
02 

(3-60) 

and 

R  '  k  n 

2  _  *2  k 1 

r2  r2p; 

(3-61) 

Clearly,  k2  cannot  vanish  if  {c..}  is  to  define  a  three 
nate  system.  The  above  equations  imply 

dimensional  coordi 

dr^ 

]  r2  R' 

h  =  [p 

cos  e]k2 

(3-62) 

or  equivalently 

/  dr  "l- 
'  r^R 

C2  =  e  1 

cos  6  . 

(3-63) 

Given  the  non-trivial  <fi  dependence  hypothesized  of  r, . , 

Theorem  2  forces 

=  r  sin  e  cos 

, 

(3-54a) 

C2  =  r  cos  6  , 

(3-64h) 

and 

;3  =  r  sin  e  sin 

4>  • 

( 3-64  c  ) 
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This  is  the  desired  conclu- 


A  similar  result  follows  if  one  takes  — 2-  =  0. 
si  on. 

Theorems  2  and  3  rather  severely  limit  the  extent  to  which  it  is 
possible  to  generalize  the  dipolar  field  aligned  coordinate  system  in  the 
manner  described  in  Section  2,  using  separable  orthogonal  coordinates.  To 
see  this,  note  that  such  a  generalization  must  involve  one  coordinate 
function  Ci(r)>  say,  which  satisfies 


^(r)  +  VdiPole  =  —‘V1  for  vaU,es  of  r  '  Rearth  (3‘55) 

-► 

Here,  p  is  the  moment  vector  associated  with  the  dipole  field  to  which  the 
coordinate  system  involving  Ci  is  aMgned  in  the  near  earth  region,  how 

■>  * 

p  .  r  =  pj  r  sin  e  cos  4>  +  p2  r  sin  e  sin  $  +  p 3  r  cos  9  (3-hh) 

Jr 

In  order  for  the  coordinate  function  d  to  be  separable,  p . =  p.  e.  for 
some  particular  i.  If  i  =  1  or  2,  then  Theorem  3  implies  that  Ci  rust  be 
part  of  a  cartesian  system  and  so  cannot  become  dipolar  for  any  value  of 
r.  If  i  =  3,  then  Theorem  2  completely  fixes  the  set  of  coordinate  func¬ 
tions  to  which  belongs  once  the  radial  dependence  of  is  specified. 

In  what  follows,  the  coordinate  system  resulting  from  the  choice 
of  Ci  given  by 

Cl  =  (r  -  kr  n)  cos  9  ( 3-£> 7 ) 

for  arbitrary  non-zero  constants  k  and  n  will  be  investigated  with  partic¬ 
ular  emphasis  upon  the  case  n  =  2.  For  large  r  and  positive  n,  Cj 
approaches  the  Cartesian  coordinate  z.  For  small  r  and  n  =  2,  C1 
approaches  the  dipole  coordinate  cos  a.  Note  that  the  spherical  surface 


of  radius  k  ,  centered  at  the  origin  is  a  surface  of  constant  Thus, 

♦ 

the  orthogonal  coordinate  system  including  Cj(r)  also  possesse:  the  quali- 
tative  features  of  a  spherical  coordinate  system  for  a  certain  range  of  r. 


f— 


SECTION  4 

A  SPECIAL  COORDINATE  SYSTEM 

The  complete  specification  of  the  set  of  orthogonal  coordinate 
functions  to  which  (Equation  3-57)  belongs  can  be  accomplished  by  using 
Theorem  2.  To  do  this,  the  following  indefinite  integral  must  be  evalu¬ 
ated  . 


V 


i 


l« '  (_J 

r2  ( 1+nk r  n  * ) 


(4-la) 


=  j  d  r  f  -  — -  + 
nr 


n  +  1 
n 


n 


+  kn 


(4-lb) 


=  in[(r  +  kn/r)^n]  +  irrelevant  constant  .  (4-lc) 

Inserting  this  expression  into  Equations  (3-3)  above,  the  following 
coordinate  functions  are  obtained: 


41  =  (r  -  kr  n)  cose  ,  ( 4-2a ) 

42  =  (r°  +  nkr_1)1//n  Sin9  >  (4-2b) 

and 

?3  =  4> 
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I  I  ■’ 


V 


i 


*! 


Specializing  to  the  case  n=2,  Equations  (4-2)  become 


Cl  =  ( r  -  k/r2)  cos  0  , 


c2  =  /( r 2  +  2k / r )  sin  8  , 


and 


43  "  ^ 


(4-3a ) 
(4-3b) 

( 4-3c ) 


As  a  check,  the  orthogonality  condition 


vCi  •  =  0  ,  i  *  j  , 


(4-4) 


may  he  explicitly  verified  hy  computing  the  gradients  of  the  4.  above: 


2k 


V4t  =  (1  +  — =■)  cos  9  r  -  (1  -  k/r3)  sin  e  9  , 


(4-Sa ) 


?42  =  ^-r  k/rJ-  sin  9  r  +  i  2(r2  +  2k/r)  cos  0  9  ,  ( 4 -b h ) 

✓r2  +  2k /r 


V43  =  - l -  l  .  ( 4- 5c  ) 

r  sin  9 


1 


While  the  orthogonality  of  the  coordinate  system  defined  hy 
{c.  }  is  guaranteed  hy  Theorem  2,  the  question  of  invertihi 1 ity  of  this 
coordinate  system,  i.e.,  of  whether  or  not  the  specification  of  (4j.C2.43) 
uniquely  determines  (r,9,<j>),  has  yet  to  be  addressed.  In  order  to  treat 
this  question,  it  is  useful  to  consider  two  separate  cases:  ci=0  and  Cj*0. 
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In  this  case,  hy  Equations  4-3,  there  are  two  alternatives.  Either: 


r  =  k1/3  , 


(4-6a ) 


6  =  sin-1  (— _ ’ll — )  ,  (4-6b) 

✓3  k1/3 

and 

4>  =  4  3  ( 4 -6  c  ) 

or 

9  =  ±  u/2  ,  ( 4 - 7 a ) 

r3  -  $2r  +  2k  =  0  ,  (4-7h) 

and 

*  =  c3  -  ( 4- 7c ) 

The  "or"  here  i^,  of  course  an  inclusive  "or",  and  the  sign  of  9  is  the 
same  as  the  sign  of  c2.  The  branch  of  this  alternative  which  obtains  for 
a  given  triple  (0,^2,^3)  (and,  hence,  the  answer  to  the  question  of  inver- 
tibility)  is  determined  by  the  magnitude  of  q2.  If  | |  >  73  k  l'  3,  then 
Equation  (4-6b)  has  no  real  solution  so  that  Equations  (4-7)  must  be 
used,  the  radial  variable  r  is  then  determined  by  the  positive  roots  of 
Equation  (4-7h).  As  discussed  in  the  Appendix,  these  roots  are  two  in 
number.  One  is  always  less  than  k1^3  and  the  other  is  always  greater  than 
the  k1^3.  They  are  given  explicitly  by10 


and 


2 

ri  =  — -  |C2lcos  A 
73 

r2  =  L  |  C2  j  cos  (A  +  y  ) 


(4-8a  ) 


(4-8b) 
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with 


.  1  _  i  r -3 «/3  k 

A  =  -  ros  - 


(4  -He  ) 


Thus,  the  5  coordinate  system  is  not  strictly  invertible  for  t1=U,  and 
| C 2 1  >  ^3  k1/3.  On  the  other  hand,  if  |  c2|  <  ^3  k  3,  then  Equation 
( 4—  7h )  has  no  real  roots  (c.f.  Equation  4-8c),  so  the  use  of  Equations 
( 4-6)  is  necessitated .  These  equations  yield  unique  values  of  (r.  e,  f>) , 
thus  insuring  invertihi 1 ity  for  Cj=0  and  |  c2j  <  *3  k  ^  3.  Finally,  if  |  t?j 
=  /3  k1/3,  then  Equations  (4-6)  and  Equations  (4-7)  give  the  sane  unique 
result,  namely  r  =  k1/3,  0  =  ±  tt/2. 


#0 


In  analyzing  this  case,  it  is  useful  to  begin  by  isolating  the 
theta  dependance  in  Equations  ( 4 -2a )  and  (4-?b): 


--1 — =  =  cos  9 
•  k/r2 


( 4- 9a ) 


- iz..—  =  sin  0  ,  ( 4 - 9h ) 

+  2k /r 

Squaning  both  sides  of  each  equation  and  adding  the  results,  there  follows 

-il  r  —  ♦  -il-  •  1  .  (4-10) 

(r 3  -  k ) 2  r3  +  2k 

This  expression  can  be  simplified  by  multiplying  both  sides  by 
(n3-k)2  (r3+2k),  expanding  the  various  products  which  arise,  and  collect¬ 
ing  powers  of  r.  Equation  (4-10)  then  becomes 
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P(r)  =  r  ^  -  Uf+^r7  -  2k(  -  3k  2r 3  -  k^r  +  2k  3 

=  0  .  (4-11) 

Equation  (4-11)  implictly  determines  r  for  given  and  42  terms  of  the 
(positive)  roots  of  a  ninth  order  polynominal  P(r).  Unfortunately,  the 
complexity  of  Equation  (4-11)  appears  to  render  dubious  the  prospect  of 
obtaining  explicit  analytic  expressions  for  these  roots  as  was  done  above 
for  the  case  ^=0.  It  is  possible,  however,  to  determine  rigorously  the 
number  of  positive  roots  of  P(r)  and  to  obtain  upper  and  lower  bounds  upon 
each  of  these  roots  without  benefit  of  such  explicit  knowledge.  This  is 
done  in  detail  in  the  Appendix,  using  certain  theorems  from  the  theory  of 
equations.  Therein,  it  is  shown  that  P(r)  has  precisely  two  positive 
roots  for  c^O.  One  of  these  lies  in  the  interval  (0,  k1^3)  and  the 
other  in  the  interval  ( k  1  / 3 ,  k1/3  +  /M )  where 

M  =nax[(c7+c \),  3k 2/ 3 ]  (4-12) 

Thus  the  4-coordinate  system  is  double  valued  for  non-vanishing  values  of 
4 1 ,  just  as  it  is  for  4^0,  42>/3  k1/3. 

Even  though  the  4-coordinate  system  is  not  strictly  invertible, 
as  shown  above,  the  lack  of  invertibi 1 i ty  is  of  a  sufficiently  innocuous 
form  that  this  coordinate  system  may  still  be  used  in  practical  applica¬ 
tion.  It  is  simply  necessary  to  specify  a  parameter  which  distinguishes 
between  the  inside  and  outside  of  the  sphere  of  radius  k1/3,  in  addition 
to  the  triple  (4x»42»C3) -  Then  the  polynomial  P(r)  can  be  numerically 
solved  for  a  unique  value  of  r,  either  lying  in  the  interval  (0,  k1,/3) 
or  in  the  interval  (k1/3,  k1/3  +  /M).  A  mesh  which  has  been  generated 
using  this  technique  is  shown  in  Figure  2.  The  value  of  k1/3  which  has 
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3i 


r.np  circular  area  in  the  middle  of  the  figure.  The  z-axi 
aligned  with  the  geomagnetic  dipole  axis  in  this  example. 


bpen  chosen  here  is  6R^ .  The  particular  nunerical  algorithm  employed  in 
solving  for  the  positive  roots  of  P(r)  was  the  Newton-Raphson  method. 

In  addition,  the  following  constraints  have  been  adhered  to  in  selecting 
the  c i  and  c2  contour  values  for  this  mesh. 

a)  Cells  on  the  surface  of  the  earth  at  a  colatitude  of  lb° 
have  dimension  -  100  km  x  100  km. 

h)  The  dimensions  of  successive  cpII  boundaries  along  the 
z-axis  moving  away  from  the  earth  increase  by  15?  until  a 
dimension  equal  to  is  reached  whereupon  the  rafp  of 
increase  is  reduced  to  zero. 

c)  In  the  upper  herd  sphere ,  the  dimensions  of  successive  cell 
boundaries  along  the  surface  of  the  earth  are  increased  by 
15';'  moving  both  clockwise  and  counter-clockwise  from  a 
rolafitude  of  15°.  boundaries  in  the  lower  hemisphere  are 
obtained  by  reflection  in  the  x-y  planp. 

d)  The  dimensions  of  successive  cell  boundaries  along  the  y- 

c.xis  to  th-  right  of  th<  circular  boundary  ^rrease  by  15° 

until  a  dimension  equal  to  R  is  attained  whereupon  the  rate 

e 

of  ircreased  is  decreased  to  zero. 

'’>■  •  1  :  r  u*  r  i  ’  r » •  ,iVe  rp*>n  imposed  in  an  attempt  to  realistically  meet 
*  •  ■  '  >  •  <  ,  *  an  Vf  -1 1  rode  des  i  lined  to  simulate  the  behavior  of  a 

■  ’  ’  i  •  r  •  cj-i>..rir  |lasma.  Such  a  code  should  he  capable  of  fine 

'*  ' 1  *  ’  ■  r  •  *  >  e  v'in'*y  rf  ft-e  auroral  oval  and  at  the  same  time  be 

■  .  *  '■i  •  1  ..  ''i-  n,  features  of  magnet  nspbp  r  i  r  current  systems.  The 

.  '  *’■  '  at  .*  i  *  Kurpiitc  rf  surrpss’ve  cell  boundaries  is 

‘‘t  re,.,-  *  rorimize  tup  total  number  of  rells  in  t be  mesh 

■'  •  •  .  » *  ,  .-.•iirt  i.omput  er  storage  required  by  the  rode)  while 


Figure  A-l.  Horner's  Method  applied  to  the  polynomal  H(R  ]  . 


The  first  substitution  transforms  the  equation  P(R)=U  into  the 
equation  P(R‘)=0  with 

P(R')  =  2R  '  9  -  2  2R  '  8  '  3R  1  6  -  R'6  -  (7  f+7  o)R  *  ?  +  1  .  (A-7) 

As  before,  implementing  the  next  two  substitutions  rears  rewriting  P(RM 
as  a  polynomial  in  R'"  =  R'-l: 

9 

P(R')  =  l  b  (Zlt  Z2)  (P‘-l)n  . 
n=n  n 

This  polynomial  may  also  be  shown  (see  b^low),  to  possess 
exactly  one  variation  of  sign  for  all  Zj  and  Z2,  Z)*C.  Again  by 

a 

Descartes'  Pule  of  Signs,  this  means  that  P(R ' )  has  one  and  only  ore  root 
greater  than  unity  for  each  Zj  and  Z2,  Zj*0.  This  in  turn  implies  that 
P ( R )  has  one  and  only  one  positive  root  less  t^an  unity.  Hence,  P(r)  has 
one  and  only  one  root  less  than  k1/3  for  each  Zj  and  Z  2,  Zj*0.  Since 
r=k'/3  is  not  a  root  of  P(r)  for  Cj*0,  the  assertion  that  P ( r T  has  only 
two  positive  roots  -  one  less  than  and  the  other  greater  than  k  3  - 
will  have  been  demonstrated  once  it  has  been  shown  that  the  coefficients 
ja^}  and  (h  },  for  0  _<  n  _<  9,  possess  one  and  only  one  variation  of  sign 
for  arbitrary  Zj  and  Z2,  Zj*0. 

One  way  to  establish  this  property  of  the  coefficients  a^,  i 
to  calculate  explicitly  these  quantities  using  Horner's  Method.  Figures 
A-l  and  A-?  show  the  results  of  applying  this  algorithm  to  the  polynom- 

~  at 

ials  P(R)  and  P(R)  respectively.  The  required  coefficients,  listed  in 
Tahlp  1  for  convenience,  are  given  hy  the  underlined  terms  in  these 
Figures.  Consider  first  the  ap's.  ag  and  ag  are  positive  while  a!  and  a 


of  thp  initial  substitutions  in  this  sequence.  One  particularly 
advantageous  choice  is  that  set  of  substitutions  which  results  in  a 
translation  of  R  by  unity.  This  set  of  substitutions  is  given  by 

R  =  —  +  1  (A -4a  ) 

R'  =  1/R "  ( A -4b ) 

which  leads  to 

R  =  R"  +  1  .  (A-4c) 

This  choice  is  intuitively  appealing  because  R=1  corresponds  to  r=k1/3, 
the  radial  value  for  which  the  character  of  the  ^-coordinate  system 
changes  from  being  dipolar  to  being  cylindrical.  More  importantly,  the 
polynomial  P ( R )  rewritten  as  a  polynomial  in  R"=R-1, 

9 

P(R)  =  l  a  (Zj,  Z  )  (R-l )n  ,  (A-b) 

„  -n  n 


may  be  shown  to  possess  exactly  one  variation  of  sign  for  all  Zj  and  Z2, 
Z^O  (see  below).  Ry  Descartes'  Rule  of  Signs,  this  means  that  P(R),  has 
one  and  only  one  root  Rq  which  satisfies  (R-l ) >0 ,  i.e.,  P(r)  has  one  and 
only  one  root  greater  than  k1/3.  To  find  the  number  of  positive  roots 
of  P(r)  less  than  k1/3,  consider  the  substitutions 


R  =  ^7  ( A -6a  ) 

R 

R’  =  ^  +  1  (A-6b ) 


R" 


(A-6c) 
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is  a  real  polynomial  with  the  first  negative  coefficient  preceded  by  k 
coefficients  which  are  positive  or  zero,  and  if  G  denotes  the  greatest  of 
the  absolute  values  of  the  negative  coefficients,  then  each  positive  root 
is  less  than  1  + 

a 


The  mathematical  armada  assembled  above  may  now  be  brought  to 
hear  upon  the  problem  of  analyzing  the  positive  roots  of  P(r).  It  is 
convenient  to  first  scale  P(r)  by  k3. 


(A-2 ) 

( A -3a ) 


and 


(A - 3h ) 


(A-3c ) 


Clearly  R  is  a  root  of  P(R)  if  and  only  if  r  e  k1/3  R  is  a  root  of 
o  oo 

P(r).  Hence,  the  objective  of  this  Appendix  can  be  accomplished  hy  study¬ 
ing  the  positive  roots  of  P(R).  As  suggested  above,  one  way  to  proceed  in 
this  study  is  to  apply  a  sequence  of  the  substitutions  prescribed  by 
Vincent's  Theorem  to  P(R)  until  one  arrives  upon  a  polynomial  with  no  more 
than  one  variation  of  sign.  The  amount  of  labor  required  to  pursue  this 
approach  can  be  minimized  by  a  judicious  {or  perhaps  fortuitous!)  choice 


48 


The  binominal  theorem  implies  that  this  polynominal  can  be  written  as  a 
polynomial  of  degree  n  in  (x-c): 

f(x)  =  An  {x-c)n  +  An-1  (x-c)"'1  +  ...  +  Ao  . 

Horner's  Method  is  an  algorithm  for  calculating  the  coefficients  A-.  It 
works  as  follows: 

First  calculate  A  .  Do  this  by  using  synthetic  division  to 
divide  (x-c)  into  f(x).  The  result  will  be  a  polynomial,  fj(x), 
of  degree  n-1  in  x  and  a  remainder  which  must  be  .  Now  divide 
fj(x)  by  (x-c)  again  using  synthetic  division.  The  result  is  a 
polynomial,  f2(x),  of  degree  n-2  in  x  and  a  remainder  which  must 
be  Aj.  This  sequence  is  repeated  for  a  total  of  n  tines.  The 
final  iteration  produces  the  obvious  result  An=an.  In  practice, 
it  is  convenient  to  arrange  the  various  polynomial  coefficients 
occurring  within  this  procedure  in  a  superdi agonal  array  .The 
ith  row  of  this  array  consists  cf  the  coefficients  of  f . (x) 
arranged  from  left  to  right  in  decreasing  rank  order  with  the 
last  term  being  the  coefficient  A..  An  example  of  such  an  array 
is  shown  in  Figure  A-l, 

The  final  result  from  the  theory  of  equations  which  will  be  required  in 
this  Appendix  is  a  theorem  which  bounds  the  positive  roots  oi  a  real 
polynomial . 

Result  4: 


If 


Hx)  =  an  xn  ♦  an_j  x"*1  ...  ♦  aQ  -  0  , 


>  0  , 
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m  is  counted  as  m  roots.  Two  consecutive  terms  of  a  real  polynomial  are 
said  to  present  a  variation  of  sign  if  their  coefficients  have  unlike 
signs . 


As  an  example,  the  polynomial  Equation  (A-la)  has  two  variations 
2  2  2  2 
of  sign  if  Ci  _>  ?2  anc1  four  variations  of  sign  if  Si  <  s2-  In  the  first 

case,  Descartes'  rule  allows  zero,  two,  or  four  positive  roots.  Clearly 

Descartes'  Rule  hy  itself  will  precisely  determine  the  number  of  positive 

roots  of  a  real  polynomial  only  if  the  polynomial  has  no  more  than  one 

variation  of  sign.  The  following  result  is  therefore  a  very  useful 

adjunct  to  Descartes'  Rule. 

Result  2:  (Vincent's  Theorem) 

If  an  equation  without  multiple  roots  is  transformed  succes¬ 
sively  hy  the  substitutions 

x  =  a+l/y  ,  y  =  b+l/z  ,  z  =  c+l/t  ... 

where  a,  b,  c,  ...  are  arbitrary  positive  integers,  the  transformed  equa¬ 
tion,  after  a  sufficiently  large  number  of  such  transformations,  possesses 
either  no  variations  of  sign  or  just  one. 

Descartes'  Rule  of  Signs  and  Vincent's  Theorem  form  the  basis  of 
a  universally  applicable  method  of  isolating  the  positive  roots  of  a  real 
polynomial.  In  using  this  method  it  is  helpful  to  have  an  algebraic 
algorithm  for  efficiently  performing  the  transformations  required  by 
Vincents'  Theorem.  One  such  algorithm  is  known  as  Horner's  Method. 

Result  3:  (Horner's  Method) 

Suppose 

, .  .  n  n-1 

f(x)  =  ap  x  +  an_j  x  ...  +  aQ  . 
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APPENDIX 

In  this  Appendix,  it  will  he  shown  that  the  polynomial  P(r) 

given  by 

P(r)  =  r9  -  -  2k(  -  3k2r3  -  k  zZ,\r  +  2k  3  (A-la) 

has  precisely  two  positive  roots  if  S^O.  One  of  these  roots  lies  in  the 
interval  (0,  k1/3)  and  the  other  in  the  interval  (k  V  3  +  /f’)  where 

M  =  roxfuJ+C*)  ’  2Uf-Cz)  ,  3k 2/3]  . 

As  an  immediate  corollary  of  this  result,  it  will  also  he  shown  that  the 
polynomial 

0 ( r )  =  r3-?2  r  +  2k  (A-lb) 

has  precisely  two  positive  roots  for  >  3k2/3.  One  of  these  is  always 
less  than  k1/3  and  the  other  always  greater  than  k  V 3.  The  demonstration 
depends  crucially  upon  several  results  from  the  theory  of  equations  which 
are  stated  below  without  proof.  Proofs  of  these  results  may  be  found  in 
references  10  and  11 . 

Result  1:  "Descartes'  Rule  of  Signs" 

The  number  of  positive  reel  roots  of  a  polynomial  with  real 
coefficients  is  either  equal  to  the  number  of  its  variations  of  sign  or  is 
less  than  that  number  by  a  positive  even  integer.  A  root  of  multiplicity 
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in  this  scheme  is  generated  using  the  ^-coordinate  system  and  therefore 
possesses  the  advantages  and  disadvantages  described  in  the  preceding 
paragraph,  i.e.,  it  is  well  tailored  to  the  elongated  geometry  of  the 
magnetosphere  and  allows  simple  solar  wind  boundary  conditions,  but  it 
also  requires  the  existence  of  a  spherical  interface  which  would  probably 
require  special  treatment  within  an  MHD  simulation  code. 


4Z 


accomodates  both  the  dipolarity  of  the  near  earth  magnetic  field  and  the 
elongated  structure  of  the  magnetosphere  at  large  distances.  An  example 
of  this  rejoining  is  shown  in  Figure  3c.  This  particular  nesh  may  he 
completed  in  three  dimensions,  under  the  previously  stated  constraints 
regarding  the  azimuthal  dimensions  of  cells,  by  setting  the  angular  separ¬ 
ation  between  successive  azimuthal  planes  of  the  inner  and  outer  meshes  to 
0  .Oh  and  0.07  radians,  respect i vely .  F'otice  that  the  use  of  an  azimuthal 
axis  aligned  along  the  earth-sun  line  allows  the  outer  mesh  to  be  much 
shorter  in  the  sunward  direction,  than  in  the  anti-sunward  direction 
correspondi ng  to  the  fact  that  the  magnetospherir  bow  shock  region  is 
compressed  toward  the  earth  while  the  magnetotail  is  distended  away  from 
the  earth.  Also,  the  use  of  such  an  axis  permits  a  simple  specification 
of  solar  wind  boundary  conditions.  An  apparent  disadvantage  of  the  use  of 
this  type  of  hybrid  mesh,  for  purposes  of  constructing  an  MHD  simulation 
code,  is  the  presence  of  a  spherical  surfacp  (at  r  =  k1^3)  upon  which  the 
boundaries  of  adjacent  cells  are  not  aligned.  This  surface  would  probably 
rpquire  special  treatment  within  such  a  code  if  spurious  numerical  effects 
are  to  be  avoided. 

There  are  other  ways  to  construct  "hybrid"  meshps  using  the 
^-coordinate  system  which  should  be  mentioned  briefly  here  for  complete¬ 
ness.  For  example,  instead  of  using  the  z’  axis  shown  in  Figure  3b  as  the 
azimuthal  axis  of  the  exterior  mesh,  one  might  use  the  z-axis  of  the 
^-coordinate  system  itself,  now  aligned  along  the  magnet,  ospheri  c  axis  of 
elongation.  This  choice  would  concentrate  resolution  in  the  equatorial 
regions  of  the  mesh  at  the  expense  of  resolution  in  the  polar  regions  so 
that  it  does  not  seem  to  be  as  desirable  an  alternative  as  that  pictured 
in  Figure  3  for  the  high  latitude  problem  of  concern  here.  Another  option 
is  to  generate  the  interior  mesh  of  Figure  3  using  a  spheri cal  coordinate 
system.  The  benefits  of  dipolarity  are  thereby  relinquished  in  favor  of 
the  simplicity  of  a  spherical  geometry  near  the  earth.  The  exterior  mesh 
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maintaining  acceptable  hounds  upon  numerical  error  arising  from  off-center 
differencing.  To  complete  the  mesh  in  three  dimensions,  the  contours  of 
of  Figure  2  must  be  rotated  azimuthally.  In  order  for  the  azimuthal 
dimensions  of  cells  at  the  periphery  of  the  mesh  to  be  limited  to  -  Rp 
and  to  maintain  100  km  resolution  in  the  vicinity  of  the  auroral  oval,  the 
angle  between  successive  azimuthal  planes  in  the  full  three  dimensional 
mesh  should  be  taken  equal  to  ~  0.05  radian. 

The  three  dimensional  mesh  described  above  fails  to  make  optimal 
use  of  the  cylindrical  character  of  the  ^-coordinate  system  at  large 
distances  because  the  z-axis  of  this  mesh  is  aligned  with  the  geomagnetic 
dipole  axis  rather  than  with  the  earth-sun  line.  This  problem  may  be 
alleviated  by  using  the  ^-coordinate  system  to  generate  a  mesh  with  two 
di fferent  azimuthal  axes.  One  way  to  do  this  is  illustrated  in  Figure  3. 
In  Figure  3a,  the  dipolar  part,  of  the  mesh  of  Figure  2,  i.e.,  that  part 
which  lies  within  the  sphere  of  radius  k1/3,  has  been  isolated  and  is 
shown  with  its  azimuthal  axis  directed  along  the  geomagnetic  dipole  axis, 
just  as  it  is  in  Figure  2.  In  Figure  3b,  on  the  other  hand,  the  remainder 
of  the  mesh  exterior  to  the  sphere  of  radius  k1/3,  is  pictured  with 
another  axis,  called  z‘,  being  used  to  define  the  azimuth,  z'  is  the 
azimuthal  axis  of  a  coordinate  system  which  is  related  to  the 

^-coordinate  system  by  a  simple  colatitudinal  translation: 

C*  =  (r  -  k/rz)  sin  9  ,  (4-13) 

1 

C '  =  (r2  +  2k/r)  cos  9  ,  (4-14) 

2 

and 

t '  =  <t>  •  (4-15) 

3 

Since  the  meshes  of  Figures  3a  and  3b  both  possess  a  spherical  boundary  at 
r  =  k 1 / 3 ,  they  may  be  easily  rejoined  to  create  a  "hybrid"  mesh  which 


Table  1.  The  coefficients  {an},{bn}  calculated  by  Horner’s  Method. 


2 

ao  ~  "  3Z  1 

b0  =  -3ZJ 

dl  =  -15Z  j 

hi  =  -m\ 

a2  =  -9(Z^Za)  -  24Z 1  +  27 

b2  =  27-21Z  1  - 

9Z  2 

a3  =  -27(Z f+Z |)  -  16ZJ  +  81 

b3  =  108  -  20Z  j 

-36Z  2 

a4  =  -33 ( Z ^+Z ^ )  -  4Z*  +  126 

b4  =  207  -  lozj 

-  60Z  2 

a  5  =  +  126 

b  5  =  234  -  2Z  j  - 

■  54Z  2 

a6  =  -7 ( Z  i+Z 2 )  +  84 

b6  =  165  -  28Z  2 

a  7  =  -(Z  f+Z  |)  +  36 

b7  =  72  -  8Z 2 

a8  =  9 

h8  =  18  -  Z 2 

a9  =  1 

bg  =  2 

are  negative  for  Zj*0.  Thus,  there  is  always  at  least  one  variation  of 
sign  presented  by  these  coefficients.  It  is  easy  to  see  that  there  are  no 
more  variations  of  sign  for  arbitrary  Zl  and  l2  with  ZjfO  by  considering 
possibilities  (a)-(f),  below,  regarding  the  magnitude  of  the  quantity 
Zi+Z2. 
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(a)  l\*l\  >  36 

In  this  case  a7  is  obviously  negative.  In  addition,  each  of 
the  terms  a  ,  n<7,  are  also  negative  so  that  there  is  only 
one  variation  of  sign  presented  by  the  coefficients  {a ^  }  for 
this  case. 

(b)  12  <  Z1+Z2  1  36 

Mow  a7  is  non-negative,  but  a^  for  n£  6  is  negative  defi¬ 
nite.  Again  there  is  only  one  variation  of  sign. 

(c)  6  <  l\+l\  <  12 

a7  and  a6  are  both  non-negative  while  a^  for  n_<5  is  negative 
definite.  There  is  only  one  variation  of  sign. 


(d)  w  K  Z^+Z2  -  6 

a7,  a6,  and  a5  are  all  non-negative  while  an  for  n<4  is 
negative  definite.  There  is  only  one  variation  of  sign. 


(p) 


3  <  Zi+Z 


2 

2 


< 


126 

33 


a7,  a6,  a5,  are  non-negative  while  a3  and  a2  are  negative 

definite.  Hence,  regardless  of  the  sign  of  a4,  there  is 

only  one  variation  of  sign  presented  by  the  sequence  of 

coefficients  (a  }. 

1  n  1 


(f)  0  <  Zj+Z2  <  3 

In  this  case,  an  for  n_>4  is  positive.  Also 

a  3-3  2  =  -IB ( Z  i+Z  2 )  +  BZ2!  4  54  >  0  , 

so  again  there  can  only  be  one  variation  of  sign. 
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2  2 

For  arbitrary  Zt  and  Z2  with  Zj*f)  thp  quantity  (Z  ,+Z  2)  must  fall 

within  one  of  the  intervals  (a)  -  (f).  Hence  it  has  heen  shown  that  the 

sequence  of  coefficients  {a^}  presents  only  one  variation  of  sign  for  all 

Zy,  Z2,  Z:*0.  How  consider  the  coefficients  {h^}.  h0  and  bj  are  negative 

while  h g  is  positive,  so  that  there  is  at  least  one  variation  of  sign 

presented  hy  the  coefficients  of  P.  To  see  that  there  is  in  fact  only 

one  variation  of  sign  presented  hy  these  coefficients,  proceed  as  above  to 

2  2 

consider  several  possibilities  regarding  the  value  of  Z2,  with  Zi  arbi¬ 
trary  hut  non-vanishing. 


(a)  Z2  >  18 

In  this  case,  b>n  for  n^8  is  negative  so  that  there  is 
clearly  only  one  variation  of  sign. 


(h)  9  <  Z2  <  18 

Nov/  hg  and  b8  are  non-negative  while  h>n  for  n<7  is  negative 
definite.  Again  there  is  only  one  variation  of  sign. 


(c) 


165 

28 


<  Z,  <  9 


hg,  b8  and  b7  are  non-negative  and  h^,  for  n£6  is  negative; 
there  is  only  one  variation  of  sign. 


(d) 


3  <  Z 


2  165 

2  -  28 


hg,  b8,  b7,  and  h6  are  non-negative  while  b^  for  n<3  is 
negative  definite.  Also 

b5-b4  =  27  +  8Z*  +  6Z2  >  0  , 

so  that  b5>b4  for  all  Zlt  Z2.  Hence,  there  can  be  only  one 
variation  of  sign  for  this  case  also. 


(e)  l\  <  3 

For  this  case,  the  b^ 's  satisfy  the  inequality 
h5  >  b4  >  h3  >  h2  . 

The  first  part  of  this  chain  has  already  heen  established 
above.  To  demonstrate  the  remaining  parts,  compute  the 
differences 

b4-b3  =  99  +  10Z 3  -  24Z2 
and 

b3-h2  =  81  +  l\  -  27 Z 2  > 

2 

which  are  both  positive  definite  for  Z2<:3  and  Zj*0.  With  this  inequality 

2 

and  the  observation  that  hg,  bg,  b?,  and  bg  are  positive  for  Z ^  3,  it  is 

evident  that  there  is  only  one  variation  in  sign  in  the  coefficients 

lb  }  for  this  case  also. 

1  n  1 


It  has  now  been  shown  that  P(r)  has  precisely  two  positive  roots 
for  Zt*0.  One  of  these  roots  lies  within  the  interval  (0,  k1/3),  and  the 
other  lies  outside  of  this  interval.  An  upper  bound  for  the  larger  root 
ran  be  obtained  by  using  Result  4  above.  It  immediately  follows  by 
inspection  of  P ( R )  that  the  larger  root  must  lie  in  the  interval 
(p 1/3,  p  1/3  +  )  with 

M  s  maxU5+£2,  3k2/3)  . 

This  is  one  of  the  results  which  was  to  be  demonstrated  in  this  Appendix. 


To  derive  the  other  result,  namely  that  the  polynomial  0(r) 
(Fquation  A-lb)  has  precisely  two  postive  roots  -  one  less  than  and  the 
other  greater  than  k1/3,  if  |£2|  >  /3  k1/3,  note  first  that 


Cl=0 


(r3-k ) 2  0(r) 


5fc 


P(r) 


( A -9 ) 


This  follows  by  direct  multiplication  or  by  inspection  of  Fquatinn  { 4 ~  1 0 ) 
in  the  text.  Equation  (A-9)  can  be  rewritten 


P(r) 


1  Ci=0 


[ r-k 1  / 3 ) 2  ( r2+rk  3/3+k  2/ 3)  2  Q(  r) 


( A  - 10 ) 


Since  the  polynomial  comprising  the  middle  factor  of  the  right  hand  side 
of  Equation  (A-10)  has  no  variation  of  sign,  it  can  have  no  positive 
roots.  Therefore,  the  positive  roots  of  Q(r)  are  the  same  as  the  posi¬ 
tive  roots  of  P(r)j^  _Q  with  the  factor  (r-k1^3)2  removed.  Thus,  the 
above  assertion  regarding  0(r)  can  be  demonstrated  by  showing  that  the 
coefficients  {ap  }n_2  9  anc<  {bn  }n=2  9  ahove  Prpsent  one  and  only  one  vari¬ 
ation  of  sign  for  ^=0,  |52|>/3  k1/3.  That  this  is  so  can  be  seen  immedi¬ 
ately  from  the  above  discussion  regarding  these  coefficients. 
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